%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [images,labels] = loadMNIST(imageFile,labelFile,id)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% LABELS
% tic;
fid = fopen(labelFile,'r');
labelHeaderSize = 8;
header = fread(fid,labelHeaderSize,'uchar');
labels = fread(fid,inf,'uchar');
fclose(fid);

% SUBSET?
if (nargin==3)
  idx = find(labels==id);
  labels = labels(idx);
%  fprintf(1,'Loading %d handwritten %ds\n',length(labels),id);
else
  idx = [1:length(labels)]';
%  fprintf(1,'Loading %d handwritten digits\n',length(labels));
end;

% IMAGES
side = 28;
height = 20;
width = 20;
nPixel = side*side;
nImage = length(labels);
images = zeros(height*width,nImage);

% SLURP
colormap hot;
fid = fopen(imageFile,'r','b');
imageHeaderSize = 16;
header = fread(fid,imageHeaderSize,'uchar');
diffIdx = [idx(1)-1; diff(idx)-1];
for ii=1:length(labels)
  slurp = fread(fid,nPixel*diffIdx(ii),'uchar');
  img1 = fread(fid,nPixel,'uchar');
  img1 = reshape(img1,side,side)'/255.0;
  img2 = deskew(img1);
  img3 = crop(img2,height,width);
  subplot(1,3,1); imagesc(img1); axis equal tight; colormap hot;
  subplot(1,3,2); imagesc(img2); axis equal tight; colormap hot;
  subplot(1,3,3); imagesc(img3); axis equal tight; colormap hot;
  drawnow; 
  pause;
  images(:,ii) = img3(:);
end;
fclose(fid);
labels = labels';
% toc;
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DESLANT

function out = deskew(in)

% GRID
[nrow ncol] = size(in);
[x,y] = meshgrid(1:ncol,nrow:-1:1);
x = x(:); 
y = y(:); 

% TRANSLATE TO ORIGIN
m = in(:)/sum(in(:));
ux = m'*x;
uy = m'*y;
x = x-ux;
y = y-uy;

% MOMENTS OF INERTIA
ixx = m'*(x.*x);
iyy = m'*(y.*y);
ixy = m'*(x.*y);

% DESLANT
dx = ux+(x+ixy*y/iyy);
dy = uy+y;
y = uy+y;
x = ux+x;

% ROUND
dx = round(dx);
dy = round(dy);

% REIMAGE
out = zeros(size(in));
ii = find(dx>=1 & dx<=ncol);
jj = sub2ind(size(out),dy(ii),dx(ii));
kk = sub2ind(size(in),y(ii),x(ii));
out(jj) = in(kk);
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CROP

function out = crop(in,height,width)

[h,w] = size(in);
top = h-nnz(cumsum(sum(in,2)));
bot = h-nnz(cumsum(flipud(sum(in,2))));
out = in;

while (h>height)
  if (top>bot)
    out = out(2:end,:);
    top = top-1;
  else
    out = out(1:end-1,:);
    bot = bot-1;
  end;
  h = h-1;
end;

left = w-nnz(cumsum(sum(out,1)));
rght = w-nnz(cumsum(fliplr(sum(out,1))));

while (w>width)
  if (left>rght)
    out = out(:,2:end);
    left = left-1;
  else
    out = out(:,1:end-1);
    rght = rght-1;
  end;
  w = w-1;
end;
  
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


